% Routine to calculate standard errors by finding numerical derivatives of moments wrt parameter vector.

function [results SE] = sterrors(params, ~, simmom, nK, nL, nA, nF, extraA, Lmode, FC, momlist, ~)     

load('empmom.mat', 'VY'); 

D         = zeros(10-sum(params([1:8,12:13])==0),length(momlist),4);
stepmom   = zeros(10-sum(params([1:8,12:13])==0),length(simmom),4);
estimuse  = zeros(10-sum(params([1:8,12:13])==0),4);

% Start loop over all non-zero elements in vector of parameter estimates.

jj = 1;
for j = [1:8,12:13]
 if params(j) ~= 0
  
  % Start by evaluating moments with a 1.0 % increase in parameter j.  
  estimuse(jj,1)      = 1.01.*params(j); 
  if estimuse(jj,3)   == params(7) && estimuse(jj,3) > 1; 
   estimuse(jj,3)     = 0.99; 
  end
  paramsuse           = params;
  paramsuse(j)        = estimuse(jj,1);
  [~, stepmom(jj,:,1)] = simulation(paramsuse, nK, nL, nA, nF, extraA, Lmode, FC, 'Search', momlist);
  
  % ... then at a 1.0 % decrease ...
  estimuse(jj,2)      = 0.99.*params(j);
  paramsuse           = params;
  paramsuse(j)        = estimuse(jj,2);
  [~, stepmom(jj,:,2)] = simulation(paramsuse, nK, nL, nA, nF, extraA, Lmode, FC, 'Search', momlist);
  
  % ... then at a 5.0 % increase ...
  estimuse(jj,3)      = 1.05.*params(j); 
  if estimuse(jj,3)   == params(7) && estimuse(jj,3) > 1; 
   estimuse(jj,3)     = 0.99; 
  end
  paramsuse           = params;
  paramsuse(j)        = estimuse(jj,3);
  [~, stepmom(jj,:,3)] = simulation(paramsuse, nK, nL, nA, nF, extraA, Lmode, FC, 'Search', momlist);
  
  % ... then at a 5.0 % decrease ...
  estimuse(jj,4)      = 0.95.*params(j);
  paramsuse           = params;
  paramsuse(j)        = estimuse(jj,4);
  [~, stepmom(jj,:,4)] = simulation(paramsuse, nK, nL, nA, nF, extraA, Lmode, FC, 'Search', momlist);
  
  % Calculate the impact of the change in parameter jj:
  [~, D(jj,:,1), ~] = evaluation(stepmom(jj,:,1),simmom,VY,momlist); D(jj,:,1) = D(jj,:,1)/(estimuse(jj,1)-params(j));
  [~, D(jj,:,2), ~] = evaluation(stepmom(jj,:,2),simmom,VY,momlist); D(jj,:,2) = D(jj,:,2)/(estimuse(jj,2)-params(j));
  [~, D(jj,:,3), ~] = evaluation(stepmom(jj,:,3),simmom,VY,momlist); D(jj,:,3) = D(jj,:,3)/(estimuse(jj,3)-params(j));
  [~, D(jj,:,4), ~] = evaluation(stepmom(jj,:,4),simmom,VY,momlist); D(jj,:,4) = D(jj,:,4)/(estimuse(jj,3)-params(j));
  [~, ~, W]         = evaluation(stepmom(jj,:,4),simmom,VY,momlist);
  
  jj = jj + 1;
  save('varcovar.mat')
  fprintf('Completed parameter vector element: %2.0f \n', j)
 end
end
% Compute the Numerical Standard Errors = (D'*Omega*D)^(-1)

V(:,:,1) = (1+1/10)*(D(:,:,1)*W^(-1)*D(:,:,1)')^(-1);
V(:,:,2) = (1+1/10)*(D(:,:,2)*W^(-1)*D(:,:,2)')^(-1);
V(:,:,3) = (1+1/10)*(D(:,:,3)*W^(-1)*D(:,:,3)')^(-1);
V(:,:,4) = (1+1/10)*(D(:,:,4)*W^(-1)*D(:,:,4)')^(-1);
nse(1,:) = sqrt(diag(V(:,:,1)));
nse(2,:) = sqrt(diag(V(:,:,2)));
nse(3,:) = sqrt(diag(V(:,:,3)));
nse(4,:) = sqrt(diag(V(:,:,4)));
NSE      = mean(nse,1);
    
ii = 1;
for i = 1 : length(params)
 if (i >= 9 && i <= 11) || params(i) == 0
  SE(i) = NaN(1,1);
  results(i,1) = {''};
  results(i,2) = {''};
 else
  SE(i) = NSE(ii);
  results(i,1) = {num2str(params(i),'%6.4f')};
  results(i,2) = {num2str(SE(i),'(%6.4f)')};
  ii    = ii + 1;
 end
end

save('varcovar.mat');

end

